Case2 - Group1

108024522 劉軒成 , 108024503 莊仕祺 , 108024466 劉倍銘 , 108024703 陳昱瑋 , 108024520 林敬皓 , 108024507 張文騰
2021/4/27

Probelm formulation

Problem 1

原始問題：能否收到物料的COA時就知道inline會異常？

資料介紹

inline為反應變數、CoA為剩下變數，表示為 \[ \{y_i, x'_i\}_{i=1}^n \] ，而\(x'_i = (x_{CoA1,...,CoAp})\)（p為變數總個數、n為sample size）。

問題公式化

考慮迴歸模型： \[ y_i = f(x'_i; \bf{\beta}) + \epsilon_i \] 可以注意到雖然把\(x'_i\)整個向量擺入模型之中，但不必要所有維度都被使用。

而我們的目標是用Maximum likelihood estimator求出\(\bf{\hat\beta}\)，進而得到 \[ \hat{y_i} = f(x'_i; \bf{\hat{\beta}}) \] 再進一步從\(\hat{y_i}\)得到 \[ \hat{z_i} = \begin{cases} 0, \text{if} \ \hat{y_i} \in R_0\\ 1, \text{if} \ \hat{y_i} \in R_1\\ \end{cases} \] \(R_0\)和\(R_1\)的選取是data driven，以下介紹：

給定一個互斥的\(R_0 \& R_1\)會決定\(S_0 \& S_1\)，分別為判定是正常的集合以及判定是異常的集合，而目標是找到至少一組\(R_{0j} \& R_{1j}\)，使得 \[ \forall i \neq j,\ (1-\frac{|S_{1i} \cap S_{N}|}{|S_{1i}|}) + (\frac{|S_{0i} \cap S_{N}|}{|S_{0i}|} + \frac{|S_{1i} \cap S_{A}|}{|S_{1i}|}) \leq (1-\frac{|S_{1j} \cap S_{N}|}{|S_{1j}|}) + (\frac{|S_{0j} \cap S_{N}|}{|S_{0j}|} + \frac{|S_{1j} \cap S_{A}|}{|S_{1j}|}) \] （\(S_N\)為真實為正常的集合、\(S_A\)為真實為異常的集合）。

Problem 2

原始問題：每天收到這麼多批原物料的COA資料，能不能從中盲測出哪一批料會有問題？

資料介紹

CoA為所有變數，表示為 \[ \{x'_i\}_{i=1}^n \] ，而\(x'_i = (x_{CoA1,...,CoAp})\)（p為變數總個數、n為sample size）。

以下有兩種將問題公式化的面向，簡言之，一為計算各點離群程度，另一為分群。

問題公式化

I:

定義一函數 \[ g:x'_i \mapsto g(x'_i) \in \mathbb{R} \] 和一閾值\(T_g\)。

假若 \[ g(x_j) \geq T \] ，則判定為異常點。需要特別提醒的是\(T_g\)的選取要使得 \[ |\{x: g(x) \geq T_g, \ x \in \{x'_i\}_{i=1}^n \}| >> |\{x: g(x) < T_g, \ x \in \{x'_i\}_{i=1}^n \}| > 0 \]

當一資料點\(x'_{new}\)新進來，即使用以上的\(g(.)\)及\(T_g\)，同樣地判斷是否 \[ g(x'_{new}) \geq T \]

II:

定義一函數 \[ h:(S_i,S_j) \mapsto h(S_i,S_j) \in \mathbb{R} \] ，\(S_i\)為一非空集合（顯然可以只有一個元素）、\(h(S_i,S_j)=0,\ \text{if} \ i = j\)。

分群的迭代步驟如下，新的群\(S_{ij} = S_i \cup S_j\)，如果 \[ \forall i,j,k,l, \ 0<h(S_i,S_j) \leq h(S_k,S_l) \] 分群直到 \[ \forall 1,...,m, \ |S_m|\geq|S_{m-1}|\geq...\geq|S_1|>1 \] 並定義異常集合\(S_A\)為\(S_A=\bigcup_{j=1}^k S_j\)和正常集合為\(S_N\)為\(S_N=\bigcup_{j=k+1}^l S_j\)，使得 \[ \frac{S_A}{S_A + S_N} \in (0,15\%] \] 該分數在這個限制下，期望愈大愈好。

當一資料點\(x'_{new}\)新進來，將它分進\(S_N\)，倘若歐式距離它最近的\(s\)個點（設定為一固定之奇數，建議為3個），命名該集合為\(E_{x_{new}}\)， \[ |E_{x_{new}} \cap S_N| \geq \frac{s+1}{2} \] 若非，則\(x'_{new}\)歸類為異常點。

Preliminary Data Analysis

Inline1

1.CoA vs Inline scatterplot

2.CoA boxplot

3.CoA vs CoA scatterplot

Inline2

1.CoA vs Inline scatterplot

2.CoA boxplot

3.CoA_8 vs CoA scatterplot

Material1

1.批次 vs CoA scatterplot

2.批次 vs missing CoA scatterplot

3.CoA vs CoA scatterplot

4.Material1 abnormal barplot

Material2

1.批次 vs CoA scatterplot

2.CoA vs CoA scatterplot

3.Material2 abnormal barplot

Material3

1.批次 vs CoA scatterplot

2.批次 vs missing CoA scatterplot

3.CoA vs CoA scatterplot

4.Material3 abnormal barplot

Material4

1.批次 vs CoA scatterplot

2.Material4 abnormal barplot

Material5

1.批次 vs CoA scatterplot

2.CoA vs CoA scatterplot

3.Material5 abnormal barplot

Missing Value Treatment

首先，我們先了解欲分析的 \(7\) 個資料集中有缺失值的CoA是那些，以及單個CoA缺失值的 總數佔總批次的比例，我們可以發現到：

  1. Inline1: 沒有缺失值
  2. Inline2: 沒有缺失值
  3. Material1: CoA_3, CoA_7 有缺失值，且都超過 \(50\%\)
  4. Material2: 沒有缺失值
  5. Material3: CoA_4, CoA_23, CoA_34 有缺失值，且CoA_4, CoA_34超過 \(50\%\)，CoA_23 僅 \(5\%\) 左右，故我們將對其補值
  6. Material4: CoA_23 有缺失值，且超過 \(50\%\)
  7. Material5: 沒有缺失值

以下我們將對Material3做補值，我們採用knn以及random forest的作法，將除了CoA_4, CoA_23, CoA_34的欄位當作input變數對CoA_23做補值

補值的結果如下圖所示，我們可以看到右邊的圖中紅色線的kNN似乎改變了CoA_23原始的分配，而藍色線的random forest則保留了原始的分配形狀，故我們將採用random forest對Material3的CoA_23做補值

Anomaly Detection

Problem 1

Linear Model

因為sample size小於變數個數，在配適線性模型之前，需要先挑選變數。此外，若考慮交乘效應，則主效應和交乘效應的變數量會遠大於sample size，因此，交乘效應的考慮將在先挑選完主效應變數後。

挑選解釋變數

採用Orthogonal Greedy Algorithm(OGA)及Forward selection，挑選至新變數使得模型中有解釋變數開始不顯著。 OGA，又名正交貪婪演算法，為一個遞迴方法，逐次挑選重要變數。

這邊是採用銀老師的方法，在OGA後加入HDIC和Trim的機制，又稱OGA三階段選模。 第一階段為OGA挑選重要變數，初次疊代，對於所有參數空間中的變數，找一個可以使log-likelihood function的絕對值一次微分可以達到最大的變數，然後做MLE估計其\(\hat\beta\)。 而後經過m\(^*\)次疊代後(註\(^*\):疊代次數在論文中有推導)，會找到m個變數和其估計之\(\hat\beta\)。

第二階段為HDIC，作為此一遞迴算則的停止法則。在第一階段中，每次疊代加入一個新變數，都會計算整個模型的loss值，則每加入變數後， loss值理想會隨之下降，直到模型加入新變數後，loss並沒有因此下降，我們就做一個HDIC的機制，讓模型選取的變數到此就停。

第三階段為Trim，進一步以HDIC為工具，對算則停止前選入的變數作更細緻的篩選。 簡單來說把第二階段找到的變數模型與loss值，當作一個baseline，而後每次從已選取的變數中，移除一個變數，與baseline做比較，因為loss是越小越好，所以當移除變數的模型的loss還比baseline還大時，代表移除的變數對模型是有影響的，所以該選入。 最後我們可以得到OGA三階段選到的重要變數。

以下為此方法的論文連結：

Orthogonal Greedy Algorithm

模型配適與診斷

Inline1

由上述介紹挑選出的模型為： \[ \hat{y}_{inline1} = 47.78 - 3.11 x_{CoA3} - 7.98 x_{CoA1} - 1.18 x_{CoA7} \]

模型的參數估計如下：

殘差符合常態檢定、配適值和殘差也沒有明顯pattern，所以認定該模型合適。

Inline2

由上述介紹挑選出的模型為： \[ \hat{y}_{inline2} = -39.56 - 0.22 x_{CoA8} - 146.15 x_{CoA13} \]

模型的參數估計如下：

殘差符合常態檢定，然而，配適值小於15之下（為正常的response值域），殘差和配適值有負斜率的線性趨勢；同樣地，配適值大於20之下（為異常的response值域），殘差和配適值也有負斜率的線性趨勢，這是因為在原始資料下，解數變數落於一定區間之內，異常點的response有相同的值，正常點也雷同，這使得該模型仍有區分正常和異常的能力，但預測值卻會失準。有鑑於配適模型的目的已達成，此模型是可以使用的。

調整判斷異常的標準

這裡提供調整判斷異常閾值的方法，此方法為data driven，並且同時使得False discovery rate (FDR)最低及Accuracy最高。該方法使用Cross validation (CV)，作法如下：在解釋變數固定下、給定某個閾值，使用leave one out進行配適，並記下落下那點在該閾值下，預測結果正確與否，依序做完所有的點，算出Accuracy和1-FDR，再換下一個閾值重作。

再以inline1的配適值與預測區間為例。可以看到單以預測值以及期望偵測任何異常的角度而言，原先給的正常區間似乎略顯寬大（唯一落在區間外的預測值是異常的點）。

inline1

透過上述方法，將閾值定在±0.31~0.39 (藍色虛線之間) ，皆能完全區分異常批號。

inline2

同樣透過上述方法，將閾值定在15~21 (藍色虛線之間) ，皆能完全區分異常批號。兩者都是因為此資料的正異常值response差異較大、sample size也較小，所以不僅閾值能容忍的區域寬大、區域內的Accuracy和1-FDR也都能各自達到1。但此現象將隨sample size變大、正異常值response靠得更近而縮小。

Problem 2

以下我們將提供三種不同的方法:

  1. Isolation Forest
  2. Hierarchical Clustering
  3. Multivariate SPC

Isolation Forest

Isolation Forest是處理Anomaly Detection問題中一種經典的算法，此演算法的中心思想為對於越異常批次的CoA越容易被決策樹用二分法的方式分開，因此當我們採用此方法時，越先被分開的批次越可能是異常的批次。以下為詳細說明：

方法假設 (Assumption)

異常資料是少數且特別的，因此符合Problem 2中對於異常CoA的定義

演算法作法 (Algorithm)

每次做二分法都是隨機從 \(p\) 個維度中選取 \(1\) 個維度，並從那個維度的全距中隨機選取一個值將資料切成兩個部分(左節點與右節點)。直到一個節點只剩一個批次後，我們可以去計算他的路徑長度 \(h(x)\)，路徑長度越短代表此資料點(批次)越為異常

以下為釋例：

從以上圖我們可以看到，\(x_o\) 比 \(x_i\) 更為特別，因此itree也用更少的直線(切分)就將 \(x_o\) 分到單一個只包含他的區域(節點)，分別對應的路徑長度為 \(h(x_0) = 4\)、\(h(x_1) = 11\)

路徑長度 \(h(x)\)

我們這次再以itree的圖像來解釋：

以上為一棵itree的示意圖，\(a, b, c, d\)為資料點(批次)，我們可以看到這棵樹總共做了 \(3\) 次切分，而 \(d\) 的路徑長為 \(1\) 最為異常，\(b,c\)為最晚切分出去，因此可能是比較正常的資料點(批次)

建造森林 (Bagging)

當我們了解如何建構一棵樹後，我們可以決定總共要種幾棵樹 \((n)\)，且每棵樹要涵蓋多少的樣本(論文建議 \(256\) 個取後放回樣本)，並透過平均路徑長的方式來決定每個資料點(批次)的異常程度

異常分數 (Anomaly Score)

最後我們將平均路徑長轉換成異常分數以將其壓縮至 \(0,1\) 之間，其轉換的方式是透過以下公式：

\[ s(x,n) = 2^{-E(h(x))/c(n)} \] 其中 \(c(n)\) 是給定 \(n\) 下的平均 \(h(x)\)

決策法則

當 \(s(x,n)\) 越接近 \(1\) 越可能為異常的資料點(批次)，我們可以主觀決定一個cutpoint \(k\)，只要 \(s(x,n) \geq k\) 即判斷該 \(x\) 為異常點

Problem 2 - Iforest Solution

我們將cutpoint \(k\) 設定為 \(85\) 分位點的 \(s(x,n)\) 且令其為 \(s^*\)，只要 \(s(x,n) \geq s^*\)，即判定該 \(x\) 為異常的批次，反之則為正常的批次

下圖 \(x\) 軸為異常分數，\(y\) 軸為批次個數，紅色線為 \(s^*\)，紅色線右方的批次即為iforest判斷的異常批次

參數設置

iforest演算法有幾個主要參數需要設置，分別為

  1. 總共需要幾棵樹 \((n)\)
  2. 建構一棵樹所需的樣本個數 (sample_size)
  3. 樹的深度 (max_depth)
  4. 隨機種子 (為了讓我們可以重新產生一樣的結果)

其中樣本數只要 \(>= 256\) 即設置樣本數 / 樹 (sample_size) \(=256\)，Iforest作者認為更大的sample_size並無法提升模型準度且耗時；種樹量 \((n)\) 均設定為 \(10000\) 以確保結果收斂；樹的深度(max_depth)設定為 \(log_2(\text{sample_size})\)， 因為此為對平均樹高的一個估計，而若我們可以在平均樹高前將某些批次分至外部節點，則這些批次越可能為異常批次，再將樹分下去也沒有意義(剩下的可能都是正常的批次)。

異常批次列表

[[1]]
[1]  1  3  4 16

[[2]]
 [1]   2   3   4   5   6   7   8  10  11  13  21  25  27  29  55  90
[17]  94  95 134 137 146 151 152

[[3]]
 [1]   1   3   4   8  10  14  16  18  19  25  29  34  36  38  39  40
[17]  41  42  48  55  63  98 102 117 124 125 126 137 184 194 215 229
[33] 233 254 256 257 268 270 271 272 273 274

[[4]]
  [1]    1    2    3    4    5    6    7    8    9   10   11   12   13
 [14]   14   15   16   19   20   23   24   26   28   34   37   44   53
 [27]   54   55   56   57   58   59   60   61   62   63   64   70   71
 [40]   72   73   75   77   78   84   86   89   91   95   99  101  102
 [53]  103  105  106  109  111  113  115  116  117  126  130  133  135
 [66]  137  138  139  140  142  143  146  150  176  178  180  201  205
 [79]  213  223  225  233  238  266  267  308  319  322  330  346  353
 [92]  361  366  370  371  373  379  381  382  384  385  386  391  393
[105]  394  401  414  425  429  430  479  484  507  525  526  527  528
[118]  531  533  534  535  541  551  554  576  614  667  683  707  708
[131]  709  710  727  731  732  754  755  756  757  758  759  760  812
[144]  813  838  929 1002 1003 1004 1005 1006 1007 1008 1009 1010

[[5]]
 [1]   1   2   3   6   9  10  11  12  19  24  41  52  62  71  82  83
[17]  94 103 111 112 113 114 115 116 121 180 185 188 204 205 207 211
[33] 212 224 239 245 246 248 252 263 264

方法驗證

此部分我們嘗試將iforest應用到Problem 1 - inline 1 & 2的資料中以確認方法的可行性：

以下將異常分數由高到低做排列，我們可以發現到inline 1的批次2為真異常且被iforest排在第二名，而inline 2的10, 11, 12為真異常，分別被排在1,7,2的位置，而這結果也會因為隨機種子的不同會得出蠻不同的結果，因為資料量過小的緣故，正常批次的群聚現象較不明顯。此比較僅供參考，但要注意iforest這方法會將CoA有出現極端值的批次視為異常批次，而此假設未必為真。

Inline 1 -

 [1]  1  2  9  8  6 10  5  3  4  7
 [1] 0.6404781 0.6177277 0.6049105 0.5904785 0.5873677 0.5660562
 [7] 0.5658629 0.5600170 0.5581988 0.5479639

Inline 2 -

 [1] 10 12  3  1  4  2 11  8  6  5  9  7
 [1] 0.6307480 0.6137182 0.6059504 0.6057707 0.6027239 0.5995936
 [7] 0.5916242 0.5907475 0.5900081 0.5857071 0.5767142 0.5617631

參考資料

以下為此方法的論文連結：

Isolation Forest Paper

Hierarchical Clustering

參數選擇

以下是分群後的結果:

使用Hierarachical 對problem 1的Inline1和Inline2來做分類，以下是分類後的結果

Inline1和Inline2兩者計算點和點之間的距離算法分別是使用歐式距離和最大距離，最大距離也就是maximum norm，是透過以下算法

\[ ||x||_{\infty} = \max\{|x_1|,\ldots,|x_n|\}. \]

用以上方法發現在Inline1和Inline2可以把異常批號和正常批號分成兩群。不過這是在知道那些是異常批號下在去計算點和點之間的距離，如果在未知的情況下，只用同一種距離方式來做計算的話是沒辦法完美的把異常批號分開的。

所以不同距離的設定對不同的Material會有不一樣的分類結果，未來在對新資料做分類時，距離的挑選也是一個必須要考慮的分析，不然可能會挑選出完全不同的批號。

異常批次列表

[[1]]
[1]  3  4 16

[[2]]
 [1]  2  3  5  6  7  8 11 13 21 25 29 55 94

[[3]]
 [1]   1   9  10  12  18  28  30  31  36  40  48  53  63  94  98 110
[17] 117 123 143 187 188 189 190 209 218 223 227 228 241 244 247 248
[33] 254 256 257 258 270 271 272 273 274

[[4]]
 [1]    3    4    7    8    9   11   12   14   15   16   19   20   24
[14]   34   37   44   54   55   56   57   58   59   60   61   62   70
[27]   71   72   73   75   77   78   84   86   89   91   95   99  102
[40]  105  106  109  111  113  115  116  117  130  133  135  137  139
[53]  140  142  143  146  176  205  213  223  225  233  238  266  267
[66]  308  319  322  330  346  353  361  381  384  484  667  707  812
[79]  813  838  929 1002 1003 1004 1005 1006 1007 1008 1009 1010

[[5]]
 [1]   1   2   3   9  10  11  12  18  20  24  41  78 112 113 116 121
[17] 165 167 180 181 185 188 195 205 207 211 212 213 224 233 239 245
[33] 246 248 252 260 263 264 265 266 268



Multivariate SPC

方法介紹

我們針對Problem 2所使用的第三種模型是多維度的SPC方法，這邊採用的是\(Hotelling's\) \(T^2\) \(shewhart\) \(chart\)。在這個control chart中，對於每一個Material，我們能透過以下定義的公式來對其中第i個批號計算統計量：
\[T_i^2=(X_i-\bar{X})^T(S^2)^{-1}(X_i-\bar{X})\]
其中\(\bar{X}\)代表該Material裡所有放入模型中的CoA的平均值，\(S^2\)則是相對應的共變異矩陣。而若以\(X_i\)代表第i個批號的所有CoA數值，我們就能透過此公式求出相對應的\(T_i^2\)值。直觀上來看，\(T_I^2\)值能夠代表第i個批號與平均值之間的差異大小，因此若某些批號的\(T^2\)值越大，就代表它越有可能屬於異常批號(OC)。這邊我們透過下方的公式來計算control limit:
\[control\ \ limit=\frac{(M-1)^2}{M}Beta_{1-\alpha}(\frac{p}{2},\frac{M-p-1}{2})\]
其中M代表此material中所含有的批號總數，p代表計算中所使用的CoA總數。當\(T_i^2>control \ limit\)發生時，我們便定義批號i為異常批號。
而在MSPC的實行過程中，我們主要考慮了兩種做法：

顯然方法1是屬於較為簡單直覺的做法，而方法2則是針對方法1做出了些微的調整。考慮兩種做法的原因主要有兩個:

考慮到上述的兩種問題，我們將會先把這兩種方法同時應用在五組Material上並進行比較，再根據結果來決定最終要採用哪種方法進行後續的分析。

方法結果比較

Material 1 Material 2 Material 3 Material 4 Material 5
方法1 16.67% 11.18% 07.58% 09.41% 07.01%
方法2 20.83% 09.87% 17.68% 13.77% 09.23%
交集 16.67% 09.21% 06.86% 09.41% 07.01%

上方表格中的數字即代表該方法在對應Material中所算出的異常批號比例，最下方則是兩種方法交集後的比例。可以發現方法2由於會將常數型CoA單獨拉出做單變量的異常批號判斷，因此使得這個方法在異常批號的辨別上較為寬鬆，也導致在大多Material當中，方法2都會傾向於得到更多的異常批號比例，甚至在Material 1中超過了20%。另外也可以發現，在兩者的交集比例中，大多的數值都會非常接近方法1所求出的異常比例，這代表方法1與方法2所找出的異常批次重複程度相當大。
接著我們將方法2中，兩個階段分別算出的異常批號比例也列進表格中：

Material 1 Material 2 Material 3 Material 4 Material 5
方法1 16.67% 11.18% 07.58% 09.41% 07.01%
方法2(階段1) 16.67% 03.95% 08.66% 08.34% 05.90%
方法2(階段2) 08.33% 07.89% 10.47% 07.85% 03.69%
交集 (階段1) 12.50% 03.29% 06.86% 03.98% 05.16%
交集 (階段2) 08.33% 07.89% 01.10% 07.85% 02.21%

方法2的階段1比例，即代表在方法2中使用常態類型的CoA，套入MSPC方法後所求出的異常批號比例；階段2的數值，則代表常數型CoA在階段2所被辨別出的異常批號比例，這兩者的交集比例就會相當於是上個表格中，方法2的異常批號總比例。而我們在這個表格中主要想觀察的，是最下方「交集」部分的階段1、2比例與方法2的階段1、2比例是否接近。在方法介紹中曾經提及，我們擔心在方法1當中，MSPC的分析結果會被少數的常數型CoA給主宰，而如果發生了上述的這個情形，則我們可以預期到表格下方的「交集」部分中，階段1的數值將會接近0%，而階段2的數值將會與方法2的階段2數值接近。換言之，這就代表我們用方法1找出的交集比例幾乎都集中在常數型CoA的異常批號上，而無法找出在其他CoA中有異常的批號，也就是被常數型CoA主宰了MSPC的表現。
但透過觀察上方表格的數值，我們可以發現前段所述的情形並未發生，並且交集後的階段1數值通常與方法2的階段1數值接近。這代表了即使我們將常數型CoA一同放入MSPC模型內，也不會發生事前所擔心的少數CoA主宰模型表現的情形。因此透過上方的方法差異分析，以及考量到找出的異常批號比例關係後，我們決定在後續分析中都只採用方法1來進行討論，也就是只考慮把所有CoA都放入MSPC模型的方法。

異常批號辨別

在決定將所有CoA都放入MSPC模型之後，我們就能透過在方法介紹中所提及的\(T^2\)以及control limit的算法，來找出5個Material各自的異常批號。詳細的結果如下方圖片所示：

圖中的橫軸即為批號的順序，縱軸為\(T^2\)的數值，因此每個黑點就代表該批號以及所對應到的\(T^2\)值。紅色虛線為算出的control limit(藍色虛線在下一小節會說明)，因此只要有黑點落在紅線上方，就代表該批號在MSPC方法中將會被定義為異常批號。而異常批號在對應Material所佔的比例如下表：

Material 1 Material 2 Material 3 Material 4 Material 5
比例 16.67% 11.18% 07.58% 09.41% 07.01%

異常批次列表

此處為MSPC方法所找出的所有異常批號列表，上至下依序為五個Material與各自所對應的異常批號。

[[1]]
[1]  1  3  4 15

[[2]]
 [1]  2  3  5  6  7  8 10 11 13 21 25 29 55 71 72 75 94

[[3]]
 [1]   1   4  12  14  17  29  31  39  40  98 102 187 188 195 222 241
[17] 244 248 268 270 271

[[4]]
 [1]    1    2    3    4    6    7    8    9   11   12   13   14   15
[14]   16   19   20   23   24   34   37   44   53   54   55   56   57
[27]   58   59   60   61   62   63   70   71   72   73   75   77   78
[40]   84   86   89   91   95   99  102  105  106  109  111  113  115
[53]  116  117  130  133  135  137  139  140  142  143  146  176  205
[66]  213  223  225  233  238  266  267  308  319  322  330  346  353
[79]  361  381  384  484  667  707  812  813  838  929 1002 1003 1004
[92] 1005 1006 1007 1008 1009 1010

[[5]]
 [1]   1   2   3   9  10  11  12  41 113 116 121 180 188 207 212 224
[17] 246 263 264

新批號預測

接著是需要判定新資料是否為異常批號的問題。在判斷現有資料是否異常時是屬於phase 1的情況，而預測新資料則是phase 2的問題。而在phase 2中，我們已經將所有異常批次移除出MSPC的模型之外了，接著只需要針對phase 1的公式做出部分的修改：

假設有新資料\(X_{i+1}\)需要進行判斷，則我們可以定義新資料所對應到的\(T^2\)值可以用下方公式計算:

\[T_{i+1}^2=(X_{i+1}-\bar{X_I})^T(S_I^2)^{-1}(X_{i+1}-\bar{X_I})\]
其中\(\bar{X_I}\)代表該Material裡所有正常批號(IC)資料的CoA平均值，\(S_I^2\)則是相對應的IC資料共變異矩陣。以\(X_{i+1}\)代表新資料的所有CoA數值時，我們就能透過此公式求出相對應的\(T_{i+1}^2\)值。

而除了平均值與共變異矩陣外，在phase 2中的control limit也有所更動:
\[control\ \ limit=\frac{p(M-1)(M+1)}{(M-p)M}F_{1-\alpha}(p,M-p)\]
其中M代表此material中所剩餘的IC批號總數，p代表計算中所使用的CoA總數。當\(T_{i+1}^2>control \ limit\)時，我們便會預測這個新批號的資料為異常批號。
需要特別注意這邊所使用的平均值、共變異矩陣與資料個數、control limit都是在移除phase 1的異常批號之後重新計算得出的結果，因此和先前方法介紹中所使用的數值皆不相同。而在預測五個Material的新資料時所使用的control limit，我們將其用藍色虛線標示在先前的五張圖片中。可以看到在phase 2中的control limit大致會與phase 1的control limit接近。而在Material 1中，紅色與藍色虛線的距離相差較遠，這主要是由於IC資料的數量過少(Material總計只有24筆數據)，因此計算出的phase 2的control limit就會與phase 1的結果差異較大。這個情況在其餘四個Material中便不會發生。

方法測試:以inline 1、inline 2為例

最後我們將上述的MSPC方法套入至inline1、inline2的數據當中做測試，觀察是否能夠正確的找出異常批號。但需要特別注意，由於在這兩筆資料中，我們的CoA個數都超過現有的批號數，因此會導致公式中共變異矩陣的反矩陣無法計算。所以這邊我們放入的CoA必須要先做篩選，而這邊我選擇的篩選方式是透過Problem 1中所提及的OGA方法，將變數從對模型影響最大至小來依序坐挑選，最後inline 1中共有9個CoA，分別是CoA3、1、7、15、4、9、14、13、6；而inline 2中共有7個CoA，分別是CoA8、13、3、10、23、15、7。

透過上述方法篩選CoA後，所畫出的批號與\(T^2\)圖形、control limit如下：

可以看到在兩組inline數據中，我們都並沒有辦法挑選出真正inline有異常的批號。而如果將inline 1、inline 2的分析結果與原始EDA的結果對照，也會發現這些被MSPC錯誤挑選出的正常批號，事實上在某些CoA中也會是異常值。而在整體資料量不足的情況下，這些批號就很容易在MSPC的方法中被挑選成為異常批號，這可能就是導致MSPC在這兩組資料中表現不佳的原因。除此之外，可能造成這個結果的原因也包含：OGA方法篩選出的CoA不恰當、有某些其他影響inline結果的CoA沒有列進資料中…等等。這部分可能就需要再針對這樣的情況進行調整。

Information Summary

Problem 1

Problem 2

各方法異常批次比較

Material 1

Material 2

Material 3

Material 4

Material 5

在剛剛的summary中主要提及了三個方法之間的異常批號關係，而接下來主要想探討交集的異常批號與它們的CoA之間的關聯。以下圖形皆代表在各個Material中，三個方法的共同異常批號裡面，有多少比例的批號在特定CoA中就已經是異常值了。這個比例值如果越大，則可以說明該CoA對於此Material的分析影響越劇烈。而直觀來看，這個比例同時也代表了若我們不進行任何分析，直接從單一的CoA來進行單變量的異常批號辨別時，有多少比例的異常批號仍然能夠被診斷出來。

CoA Importance

Material 1

在Material 1中，可以發現所有CoA的比例皆為50%，這是因為Material 1中總計只有2筆的共同異常批號，因此不容易進行判斷。

Material 2

在Material 2中，CoA6的比例達到92.3%，說明這個CoA對於我們的分析結果影響劇烈。

Material 3

在Material 3中，整體的異常批號為5筆，因此大多CoA中僅有1至2筆的差異，所以我們認為無法進行判斷。

Material 4

在Material 4中，CoA13達到62.2%。雖然不像Material 2中的CoA6能夠達到90%，但仍然能代表這個CoA對我們Material 4的分析結果有很大的影響。

Material 5

在Material 5中，CoA2、3都超過30%，CoA5則達到47.4%，都明顯高於剩下的CoA1，這代表CoA1的數值對我們Materail 5的分析影響較小。

而從這段summary來看，除了我們能夠找到對各個Material分析影響較大的CoA外，同時我們也能注意到，通常能夠得到較高比例數值的CoA，多半都是屬於常數型CoA，也就是大多批號的數值都落在某些固定值上，僅在少部分批號中有所跳動。例如剛剛提到數值特別高的Material 2的CoA6、Material 4的CoA13都屬於此類型。這代表在不同的分析方法中，這種類型的CoA都對我們的結果影響很大，因此這可能是一種需要特別留意的CoA特性。

各方法異常批次趨勢

Decision Making

Problem 1

Problem 2

批號交集表